function[x,output]=conssolve(x0,phinew,c,Y);
global P beta n m S state agent gammagrid   lowerbound rho table vs
output=[];
objval = [];
options =optimset('LargeScale','off','Display','off','MaxFunEvals',1000,'TolFun',0.00001);
[x,fval,exitflag]=fsolve(@objfun,x0,options);
output(6)=exitflag;
function f = objfun(x);
index(3)=vs;
for k=2:vs
    f(k)=x(k)/(Y(1)-sum(x(2:vs)))-((1+phinew(k))/(1+phinew(1)))^(1/rho)*c(k)/c(1);
end
f(1)=Y(1)-sum(x);
end
end
